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Symbolic algorithms for the Painleve test, special solutions, 
and recursion operators for nonlinear PDEs. 

Douglas Baldwin, Willy Hereman, and Jack Sayers 

Abstract. This paper discusses the algorithms and implementations of three 
Mathematica packages for the study of integrability and the computation of 
closed-form solutions of nonlinear polynomial PDEs. 

The first package, PainleveTest .m, symbolically performs the Painleve 
integrability test. The second package, PDESpeclalSolutions .m, computes ex- 
act solutions expressible in hyperbolic or elliptic functions. The third package, 
PDERecursionQperator .m, generates and tests recursion operators. 



1. Introduction 

The investigation of complete integrability of nonlinear partial differential equa- 
tions (PDEs) is a nontrivial matter jHGCM98 . Likewise, finding the explicit 
form of solitary wave and soliton solutions requires tedious, unwieldy computations 
which are best performed using computer algebra systems. For example, the sym- 
bolic computation of solitons with Hirota's direct method and the homogenization 
method, are covered in HN97, HZ97j- 

Recently, progress has been made using Mathematica and Maple in applying the 
inverse scattering transform (1ST) method to compute solitons for the Camassa- 
Holm equation |J03| . Before applying the 1ST method (a nontrivial exercise in 
analysis!), one would like to know if the PDE is completely integrable or what 
elementary travelling wave solutions exist. 

This is where the symbolic algorithms and packages presented in this paper 
come into play. 

In this paper, we introduce three algorithms and related Mathematica packages 
|BH03| which may greatly aid the investigation of integrability and the search for 
exact solutions. 

In Section 121 we present the algorithm for the well-known Painleve integrabil- 
ity test |AC91L ICOOL [C99 . which was recently implemented as PainleveTest .m. 
Section outlines the algorithm behind PDESpeclalSolutions .m, which allows 
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one to automatically compute exact solutions expressible in hyperbolic or ellip- 
tic functions; full details are presented in |BGH 03 . In Section 01 we give an 
algorithm for computing and testing recursion operators IHG99I ; the package 
PDERecursionOperator .m automates the steps. The latter package builds on the 
code InvariantsSymmetries .m GH97b , which computes conserved densities, 
fluxes, and generalized symmetries. 

In this paper we consider systems of M polynomial differential equations, 

(1.1) A(u(x),u'(x),u"(x),... ,u('")(x))=0, 

where the dependent variable u has M components Ui, the independent variable x 
has N components Xj , and u*^™^ (x) denotes the collection of mixed derivative terms 
of order m. We assume that any arbitrary coefficients parameterizing the system 
are strictly positive and denoted as lower-case Greek letters. 

Two carefully selected examples will illustrate the algorithms. The first example 
is the Kaup-Kupershmidt (KK) equation (see e.g. |GH97aL IW98j V 

25 

(1.2) ut^bu'^Ux + —UxU2x + 'ouu:ix + u^x■ 
^]le second example is the Hirota-Satsuma system of coupled Korteweg-de Vries 
(KdV) equations (see e.g. |AC91p . 

ut = a{&uux + u^x) - '^vvx, a> 0, 
(1-3) „ 

Vt = -3uVx - V3x- 

The computations for both examples will be done by the software. 



2. The Painleve Test 

The Painleve test verifies whether a system of ODEs or PDEs satisfies the 
necessary conditions for having the Painleve property. There is some variation in 
what is meant by the Painleve property KJH9'3. As defined in ARS80 , for 
a PDE to have the Painleve property, all ODE reductions of the PDE must have 
the Painleve property. While |ARS80] requires that all movable singularities of all 
solutions are poles, the more general definition used by Painleve himself requires 
that all solutions of the ODE are single- valued around all movable singularities. A 
later version |WTG83] allows testing of PDEs directly without having to reduce 
them to ODEs. For a thorough discussion of the Painleve property, see |C00l IG99| . 

Definition 2.1. A PDE has the Painleve property if its solutions in the com- 
plex plane are single-valued in the neighborhood of all its movable singularities. 

2.1. Algorithm and implementation. Following WT G83j . we assume a 
Laurent expansion for the solution 

CO 

(2.1) Ui(x) = g"'(x) ^-Uj,fc(x)5''(x), Uifi{x)^0 and a., e Z~ , 

k=Q 

where Ui_fc(x) is an analytic function in the neighborhood of g{x). The solution 
should be single- valued in the neighborhood of the non-characteristic, movable sin- 
gular manifold (?(x), which can be viewed as the surface of the movable poles in the 
complex plane. 

The algorithm for the Painleve test is composed of the following three steps: 
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Step 1 (Determine the dominant behavior). It suffices to substitute 

(2.2) u,(x) =u,,o(x)g"'(x) 

into to determine the strictly negative integer ai and the function Ui^o{x.). In 
the resulting polynomial system, equating every two possible lowest exponents of 
(7(x) in each equation gives a linear system to determine a^. The linear system is 
then solved for ai. 

If one or more exponents ai remain undetermined, we assign a strictly negative 
integer value to the free ai so that every equation in has at least two different 
terms with equal lowest exponents. Once ai is known, we substitute l|2.2|l back into 
1)1.1(1 and solve for Uifi{x.). 

Step 2 (Determine the resonances). For each ai and Ui_o(x), we calculate the 
integers r for which ,.(x) is an arbitrary function in (|2.1(l . We substitute 

(2.3) u.(x) = M.,o(x)5"- (x) + u,,,(x)g"'+''(x) 

into (|1.1() . keeping only the most singular terms in (7(x), and require that the 
coefficients of Mi^r(x) equate to zero. This is done by computing the roots of det Q = 
0, where the AI x M matrix Q satisfies 

(2.4) Q • U,. = 0, Ur = (lil,r U2^r ■ ■ ■ UM,rY' ■ 

Step 3 (Find the constants of integration and check compatibility conditions). 
For the system to possess the Painleve property, the arbitrariness of Ui.f.(x) must 
be verified up to the highest resonance level. That is, all compatibility conditions 
must be trivially satisfied. This is done by substituting 

I'M ax 

(2.5) u,(x)=g"'(x)^M.,fc(x)/(x) 

fe=0 

into 1(1. 1(1 . where rMax is the highest positive integer resonance. For the system to 
have the Painleve property, there must be as many arbitrary constants of integra- 
tion at resonance levels as resonances at that level. Furthermore, all constants of 
integration Mi^/j(x) at non-resonance levels must be unambiguously determined. 

2.2. Examples of the Painleve test. 

Example 2.1 (Kaup-Kupershmidt). To determine the dominant behavior, we 
substitute H2.2|l into H1.2f) and pull off the exponents of (7(x) (see TableQ]). Removing 



Term 


Exponents of g(x) with duplicates removed 


Ut 


ai — 1 




ai — 5, ai — 4, ai — 3, ai — 2, ai — 1 




2ai - 3,2ai - 2,2ai - 1 




2ai - 3,2ai - 2 




3ai — 1 



Table 1. The exponents of g(x) for H1.2|l . 
duplicates and non-dominant exponents, we are left with 
(2.6) {ai - 5,2ai - 3,3ai - 1}. 

Considering all possible balances of two or more exponents leads to ai — —2. 
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Substituting u(x) = uifi{x.)g ^(x) into H1.2|l and solving for ui^o(x) gives us 

(2.7) Mi,o(x) -24g2(x) and ui,o(x) = -3gl{x). 

For the first branch, substituting u{x) = —24:g^{x.)g^^ (x) + ui^r(x)5''^^(x) into 
H1.2|l . keeping the most singular terms, and taking the coefficient of wi.r(x), gives 

(2.8) -2(r + 7)(r + l)(r - 6)(r - 10)(r - 12).g^(x) = 0. 

Hence, r = -7,-1,6,10,12. 

While we are only concerned with the positive resonances, r = — 1 is often called 
the universal resonance and corresponds to the arbitrariness of the manifold .g(x). 
The meaning of other negative integer resonances is not fully understood KJH97j. 

The constants of integration at level j are found by substituting H2.5|l into ljl.2|l . 
where rMax = 12, and pulling off the coefficients of ^-'(x). The first few are 

24322; (x), 

6gL(x) - 8g^(x)g3:^(x) 

65L(x) - 8g^{x)g2x{x)g3x{^) + 2gl{x)g4^{x) 



(2.9) 


at j 


= 1 : 


wia(x) 


(2.10) 


at j 


= 2 : 


Ul,2(x) 


(2.11) 


at j 


= 3 : 


"i,3(x) 



The compatibility conditions are satisfied at resonance levels 6, 10, and 12. The 
remaining constants of integration uij(x) are computed but not shown here. 

Likewise, in the second branch, substitute u{x) = — 3(7^(x)g~^(x)+ui^r(x).g'"^^(x) 
into (|1.2|l . and proceed as before to get r = — 1, 3, 5, 6, 7. The constants of integra- 
tion at levels j = 1,2 and 4 are again found by substituting l|2.5(l into l|1.2(l and 
pulling off the coefficients of g-'{x). This gives, 

(2.12) atj = l: Mi,i(x) = 3g2x(x), 

/r,1o^ of ■ ( \ 3g|^W "4g:r(x)fe(x) 

(2.13) atj = 2: mi,2(x) = . 2( \ ■ 

^9x IX j 

The coefficient ui_4 at level j = 4 is not shown here due to length. At the resonance 
levels r = 3, 5, 6, 7, the compatibility conditions are satisfied and (|1.2|) passes the 
Painleve test. It is well-known (see e.g. jGHQTa, W98] ) that ()1.2|) is completely 
integrable. 

Example 2.2 (Hirota-Satsuma) . The Hirota-Satsuma system illustrates the 
subtleties of determining the dominant behavior. 

As in Example 12.11 we substitute (|2.2|l into H1.3|l and pull off the exponents 
of g(x) (listed in Table [2J|. Removing non-dominant exponents and duplicates by 



Term 


Exponents of g(x) 


Ut 


ai — 1 




ai — 3, «! — 3, ai — 3, ai — 2, ai — 2, ai — 1 




2ai - 1 




2^2 - 1 


Vt 


^2 - 1 


V3x 


a2 — 3, Q!2 — 3, a2 — 3, a2 — 2, 0:2 — 2, 02 — 1 


3uVx 


ai -1- q;2 — 1 



Table 2. The exponents of g(x) for H1.3|l . 
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term, we get {ai — 3, 2ai — 1, 2a2 — 1} from Ai and {a2 — 3, ai + a2 — 1} from A2. 

Equating the possible dominant exponents from A2 gives us a2 — 3 = ai +02 — 1 
or ai = —2. Unexpectedly, ai = —2 balances two of the possible dominant terms 
in Ai, and we are free to choose a2 such that 

(2.14) 2ai - 1 < 20:2 - 1 or - 2 < 02 < 0. 

Hence, 02 = —1 or a2 = —2. 

Using the two solutions for a.;, solving for Ui^ results in 



(2.15) 



(2.16) 



ai = -2, ui,o(x) = -4.g^(x), 
= -2, U2,o(x) = ±2^/6a5l^(x), 

fai = -2, ui,o(x) = -2gl{x.), 
]a2 = —1, U2,o(x) arbitrary. 



We substitute (|2.3|) into H1.3|l while using the results for ai and Ui^o(x). For H2.15|l . 
substituting u(x) — —Agl (x)^^^ (x) + (x)g''^^ (x) and v (x) — ±2-\/6ag^ (x)g^^ (x) 
-Kt2,r(x)(7''^^(x) into (|1.3() and keeping the most singular terms gives 



^f~ir~ 4)(r2 - 5r - 18)ag3(x) ±12V6^5^(x) \ / ui,.(x) 
^' I T4(r-4)V6^5^(x) (r - 2)(r - 7)rff3(x) ^ ^ U2,r(x) 

Setting 

det Q = -a(r + 2)(r + l)(r - 3)(r - 4)(r - 6)(r - 8)g^(x) = 

yields r = -2,-1,3,4,6,8. 

As in the previous example, the constants of integration at level j are found by 
substituting (|2.5(l into (|1.3|l and pulling off the coefficients of ^-' (x). At j = 1, 
(2.17) 

-66ag3(x) Tl2V6^.g^(x) \ f ui,i(x) \ / - 120^53 (x)g2^(x) 
Tl2V6^5?(x) 65^ (x) J y U2,i(x) y T60V6^5^(x)g2x(x) 

Thus, 

(2.18) '"i,i(x) = 4.g22;(x), M2,i(x) = ±2V6a52i;(x). 



At J = 2, 
(2.19) 



. ^ 3gL(x) ~ g:,(x)(gt(x) + 4g3:,(x)) 
"i,2(x) = 



U2,2(x) = ± 



3ffKx) 

SasLW - 4a52:(x)532,(x) - (1 + 2a)5f (x)ga;(x) 



The remaining constants of integration are omitted due to length. The compatibil- 
ity conditions at r = 3 and 4 are satisfied. At r = 6 and r = 8, the compatibility 
conditions require a — ^. 

Likewise, for H2.16(l . the resonances following from the substitution of u(x) = 
-2g2(x)g-2(x) +ui,.r(x)g''~2(x) and u(x) = U2,o(x)5"^(x) +M2,,.(x)5''-1(x) into 
()1.3(l are r = —1, 0, 1, 4, 5, 6. The zero resonance explains the arbitrariness of U2.o(x). 
Similarly, we computed all constants of integration, but ran into compatibility 
conditions at r = 5 and r — 6, which require a = i. Therefore, 11. 3|) passes the 
Painleve test if a = i, a fact confirmed by other analyses of integrability |AC91| . 
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3. Travelling Wave Solutions in Hyperbolic or Elliptic Functions 

The traveling wave solutions of many nonlinear ODEs and PDEs from soliton 
theory (and elsewhere) can be expressed as polynomials of hyperbolic or ellip- 
tic functions. For instance, the bell shaped sech-solutions and kink shaped tanh- 
solutions model wave phenomena in fluid dynamics, plasmas, elastic media, elec- 
trical circuits, optical fibers, chemical reactions, bio-genetics, etc. An explanation 
is given in jHT90j . while a multitude of references to tanh-based techniques and 
applications can be found in BGH03 , MH961. 

In this section we discuss the tanh-, sech-, cn- and sn- methods as they apply 
to nonlinear polynomial systems of ODEs and PDEs in multi-dimensions. 

3.1. Algorithm and implementation. All four flavors of the algorithm 
share the same five basic steps. 

Step 1 (Transform the PDE into a nonlinear ODE) . We seek solutions in the 
traveling frame of reference, 

N 

(3.1) ^ = J2cjX,+6, 

where the components Cj of the wave vector and the phase 5 are constants. 

In the tanh method, we seek polynomial solutions expressible in hyperbolic 
tangent, T — tanh^. Based on the identity cosh^ ^ — sinh^ ^ = 1, 

(3.2) tanh' ^ = sech^ C = 1 - tanh^ ^, 

(3.3) tanh" ^ = -2 tanh,^ + 2 tanh^ ^, etc. 

Therefore, the first and all higher-order derivatives are polynomial in T. Conse- 
quently, repeatedly applying the chain rule, 

(3.4) T = tanh(0: |l = = c,(l - T^)^ 
^ ^ dxj dT dxj ' dT 

transforms (|1.1|) into a coupled system of nonlinear ODEs. 
Similarly, using the identity tanh ^ -I- sech ^ = 1, we get 

(3.5) sech'^ = — sech^ tanh^ — — sech ^ yl — sech^. 
Likewise, for Jacobi's elliptic functions with modulus m, we use the identities 

(3.6) sn^(^; m) = 1 — cn'^(^; m), and dn^(^; m) = 1 — m + to cn^(^; m), 
to write, for example, cn' in terms of cn : 

(3.7) cn'(^;m) = — sn(^; to) dn(^; m) = — \/(l — cn2(^; to))(1 — to-|-to cn^{£,]m)). 
Then, repeatedly applying the chain rules, 

(3.8) 5 = sech(0: ^ = -'^.'^V^^^^' 

(3.9) CN = cn(C;TO) : ^ = -CjJ {I - CN2)(1 - m + m CN^), 

OXj 

(3.10) SN = sn(^;TO): |^ = 0,^(1 - SN2)(1 - mSN^)-^, 
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we transform Hl.lf) into a coupled system of nonlinear ODEs of the form 
(3.11) r(F, u(^^), u'(F), . . . ) + v/i?Mn(i^, u(F), u'(F), . . . ) - 0, 

where F is either T, S, CN, or SN, and R{F) is defined in Tabled 



F 


R{F) 


T 





S 


1-5^ 


CN 


(l-CN^)(l-m + mCN^) 


SN 


(l-SN^)(l-mSN^) 



Table 3. Values for R{F) in (jXTT|) . 

Step 2 (Determine the degree of the polynomial solutions). Since we seek 
polynomial solutions 

Mi 

(3.12) C7,(F) = ^a,,F^ 

i=o 

the leading exponents Mi must be determined before the aij can be computed. The 
process for determining Mi is quite similar to the one for finding ai in Sectional 

Substituting Ui{F) into (|3.11|l . the coefficients of every power of F in every 
equation must vanish. In particular, the highest degree terms must vanish. Since 
the highest degree terms only depend on F*^' in H3.12|l . it suffices to substitute 
Ui{F) = F^^' into l|TTT|l . Doing so, we get 

(3.13) P(F) + /r(F)Q(^^) = 0, 

where P and Q are polynomials in F. Equating every two possible highest exponents 
in each Pi and Qi gives a linear system to determine Mi . The linear system is then 
solved for Mi. 

If one or more exponents Mi remain undetermined, we assign a strictly positive 
integer value to the free M^, so that every equation in (|3.11() has at least two different 
terms with equal highest exponents in F. 

Step 3 (Derive the algebraic system for the coefficients a^). To generate the 
system for the unknown coefficients aij and wave parameters Cj, substitute (|3.12|l 
into (|3.11|l and set the coefficients of F^ to zero. The resulting nonlinear algebraic 
system for the unknown is parameterized by the wave parameters Cj and the 
parameters in (|l.l|l . if any. 

Step 4 (Solve the nonlinear parameterized algebraic system). The most diffi- 
cult aspect of the method is solving the nonlinear algebraic system. To solve the 
system we designed a customized, yet powerful, nonlinear solver. 

The nonlinear algebraic system is solved with the following assumptions: 

(1) all parameters (the lower case Greek letters) in (ILH are strictly positive. 
(Vanishing parameters may change the exponents Mi in Step 2) . To com- 
pute solutions corresponding to negative parameters, reverse the signs of 
the parameters in (|1.1() . 

(2) the coefficients of the highest power terms {aiMn i = 1, • • • , Mi) in (|3.12|) 
are all nonzero (for consistency with Step 2). 

(3) all Cj are nonzero (demanded by the physical nature of the solutions). 
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Step 5 (Build and test solutions). Substitute the solutions from Step 4 into 
(|3.12|l and reverse Step 1 to obtain the explicit solutions in the original variables, 
ft is prudent to test the solutions by substituting them into 

3.2. Examples of travelling wave solutions. 

Example 3.1 (Kaup-Kupershmidt). While the tanh-, sech-, cn- and sn-methods 
find solutions for H1.2(l . we demonstrate the steps of the algorithm using the tanh- 
method. After which, we summarize the results for the other methods. 

First, transform (|1.2|l into a nonlinear ODE by repeatedly apply chain rule 
(|S3|). The resulting ODE is 

(3.14) 2c2U' + ci llOU^U' + cl [25{T^ - 1)U'[2TU' + (T^ - l)U"] 
+ 10U[{6T^ - 2)U' + 6T{T^ - 1)U" + (T^ - i fu'" 
+ cl{16{15T* - 15T^ + 2)U' + AO{T^ ~ l)(6r(2T^ - 1)U" 
+ (6T4 - + l)U"' + T{T^ - l)2[/(4)) + - l)4t/(5))]] 1 = 0, 

where T = tanh(0 and U = U{T). 

Next, to compute the degree of the polynomial solution(s), substitute U{T) — 
r*^i into H3.14|l and pull off the exponents of T (see Table 0)). Remove duplicates 



Term 


Exponents of T with duplicates removed 


Ut 


Ml - 1 


•) 

u 


3Mi - 1 




2Mi + l,2Mi - l,2Mi - 3 


5uU3a; 


2Mi + l,2Mi - l,2Mi - 3 


U5x 


Ml + 3, Ml + I, Ml ~ 1, Ml - 3, Ml - 5 



Table 4. The exponents of T after substituting U{T) ^ T*^i. 

and non-dominant exponents, to get 

(3.15) {3Mi - 1, 2Mi + 1, Ml + 3}. 

Consider all possible balances of two or more exponents to find Mi — 2. 
Substitute 

(3.16) U{T)^aio + aiiT + ai2T^ 

into (|3.14() and equate the coefficients of to zero (where i = 0, 1, . . . , 5) to get 

(ai2 + 3c?)(ai2-h24c2) = 0, 
aii(5a?2 + 55ai2C? + 2Ac\) = 0, 
aii{bal^ci - lOaioCi + 25ai2Ci -f lQc\ + C2) = 0, 

(3.17) aii{al^ + 6aioai2 + &aioc\ - 4&ai2c\ - 24ct) = 0, 
Aaiiai2 -f AaiQai2 + \\a\iCi + 2AaiQai2Ci — 56a]^2Ci — 192oi2Ci = 0, 

lQaiQa\iCi -f 10a^Qai2Ci — 'iba\ic\ 

-80aioai2C? + SOa^aC? + 272ai2Ci + 2ai2C2 = 0. 
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Solve the nonlinear algebraic system with the assumption that ai2,ci, and C2 are 
all nonzero. Two solutions are obtained: 




(3-18) :]\ . and 



= 2cl an = 0, 
= -3ci, C2 = -cf , 

where ci is arbitrary. 

Substitute the solutions into and return to u(x, t) to get 

(3.19) u{x,t) = 16ci-24c?tanh^(cia;-176cfi + (5), 

(3.20) u{x, t) = 2cl ~ 'icl tanh^(cia; - c\t + 5). 
Using the sech-method, one finds 

(3.21) u{x, t) = -8c? + 24c? sech^(cix - ITGc^t + 5), 

(3.22) u{x, t) = ~cl + 3c? sech2(cia: - c^i + 5). 

Alternatively, the latter solutions can be found directly from the tanh-method so- 
lutions by using the identity tanh'^ ^ = 1 — sech^ ^. 
For the cn- and sn-methods, one gets 

(3.23) u{x,t) = 8c? [1 - 2m + 3mcn2(cia; - nQc\[m^ - m + l)t + (5;to)] , 

(3.24) u{x,t) = c? [l - 2m + 3mcn^(cia:; - Ci(m^ - m + l)t + 5]m)\ , 

(3.25) u{x,t) = 8c? [l + m - 3msn2(cia; - I76cl{m'^ -m+l)t + 6-m)] , 

(3.26) u{x, t) — Ci [l + m — 3m sn'^(cia; — c^(m^ — m + l)t + 6; m)] . 

Example 3.2 (Hirota-Satsuma) . As in the previous example, the tanh-, sech-, 
cn- and sn-methods all find solutions for H1.3|l . In this example, however, we will 
illustrate the steps using the sech-method. 

Transform (|1.3|l into a coupled system of ODEs, apply chain rule (|3.8|l and 
cancel the common SV^ — factors to get 

(3.27) 

C2U[-6aciUiU[-acl[{l-6S^)U[+3S{l-2S'^)U['+S^{l-S^')U'j"]+2ciU2U^ 0, 
C2U^ + 3ciUiU!^ + cl[{l-6S^)U^ + 3S{l-2S^)U^ + S^{l-S^)U^'] = 0. 

To find the degree of the polynomials, substitute Ui{S) = S^^\ 112(3) = 5*^= into 
(|3.27|l and first equate the highest exponents from A2 to get 

(3.28) Ml -I- M2 - 1 = 1 + M2 or Mi = 2. 

The maximal exponents coming from Ai are 2Afi — 1 (from the UiU[ term). Mi + 1 
(from [/("), and 2M2 - 1 (from U2U^). 

Since Mi = 2 balances at least two of the possible dominant exponents in Ai, 
namely 2Afi — 1 and Mi + 1, one is again left with 1 < M2 < Mi = 2, or 



(3.29) 
(3.30) 




Ui{S) — aio + ails' + auS^, 
U2(S) = 020 + a2iS, 

Ml = 2, [/i(5') = aio + aiiS' + ai2S'2, 

M2 = 2, 1/2(5') =020 +0215" + 022 S^. 
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To derive the algebraic system for , substitute H3.29|l into (|3.27(l , cancel common 
numerical factors, and organize the equations (according to complexity): 



(3.31) 



an 


a2ici 


= 0, 


aaiici(3ai2 


-cl) 


= 0, 


aai2Ci(ai2 - 


-2c?) 


= 0, 


a2ici(ai2 - 


-2c?) 


= 0, 


a2i(3aioCi + 


+ C2) 


= 0, 


daaioGiiCi — 2a2oa2iCi + aancj — 


ai2C2 


= 0, 


SaaliCi + 6aaioai2Ci — a^iCi + 4aai2C^ — 


ai2C2 


= 0. 



Similarly, after substitution of (|3.30() into (|3.27ll . one gets 

a22Ci(ai2 - 4c?) = 0, 

a2i(3aioCi + cf + C2) = 0, 

Ci(ai2a2i + 2aiia22 - 2a2ic\) = 0, 

Ci(3aaiiai2 — 021022 — aanc?) = 0, 
(3.32) / 9 9 2\ 

ci(3aa]^2 ^ '^22 ^ 6aai2Cj^) = 0, 

Gaaioaiici — 2a2oa2iCi + aanc? — aiiC2 — 0, 

3aiia2ici + 6aioa22Ci + 8a22C? + 2a22C2 = 0, 

3Q!a?;^ci + 6aa?Qai2Ci — a2iCi — 2a2oa22Ci + 4aai2C? — ai2C2 = 0. 

Since 012, 021, a, ci, and C2, are nonzero, the solution of H3.31|l is 

{aio = -(ci+C2)/(3ci), 020=0, 
aii=0, 021 = ±V(4act - 2(1 + 2a)ciC2), 

ai2 = 2c?. 

For ai2, 022, a, ci, and C2 nonzero, the solution of l|3.32|l is 

'aio = -(4c? + C2)/(3ci), 020 = ±[4ac? + (1 + 2a)c2]/(cl^/6^), 
(3.34) I an = 0, 021 = 0, 

ai2 = 4c?, 022 = T2c?V6a. 

The solutions of (|1.3|) involving sech are then 



uix, t) = — i h 2c? sech^(cia; + C2t + 5), 

3ci 



(3.35) 



and 



(3.36) 



civ 6a 

In both sets of solutions, ci, C2, a, and 5 are arbitrary. 



u(x, t) — ±1/ 4acj — 2(1 + 2q!)ciC2 sech(cia; + + 5), 



4c"^ C2 

u{x, t) = ^ h 4c? sech^(cia; + C2t + 6), 

, 4ac? + (1 + 2a)c2 „ 9 ,2, 
v{x, t) = ± i ^= — =F 2c?V6a sech^(cix + C2t + 5). 
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4. Recursion Operators 

In this section, we only consider polynomial system of evolution equations in 
(1 + 1) dimensions, 

(4.1) Ut{x,t) = F{u{x,t),Ux{x,t),U2x{x,t), . . . ,Ujnx{x,t)), 

where u has M components Ui and ~ d"^u/dx"^ . For brevity, we write F(u), 
although F typically depends on u and its ^-derivatives up to order m. If present, 
any parameters in the system arc strictly positive and denoted as lower-case Greek 
letters. 

The algorithm in Section will use the concepts of dilation invariance, den- 
sities, and symmetries. 

4.1. Scaling invariance, densities, symmetries. A PDE is dilation invari- 
ant if it is invariant under a dilation symmetry. 

Example 4.1 (Kaup-Kupershmidt). As an example, ()1.2|l is invariant under 
the dilation (scaling) symmetry 

(4.2) (t, x, u) (A^^t, X-^x, A^u), 

where A is an arbitrary parameter, leaving A^ as a common factor upon scaling. 

To find the dilation symmetry, set the weight of the x-derivative to one, w{Dx) = 
1, and require that all terms in (|4.1II have the same weight. For 1)1. 2|l . we have 

(4.3) w{u) -f w{Dt) = iw{u) + 1 = 2w{u) + 3 = 2w{u) + 3 = w{u) + 5, 

or u'(u) = 2 and w{Dt) — 5, Consequently, in (|1.2(l the sum of the weights or rank 
of each term is 7. 

A recursion operator^ TZ, is a linear integro-diffcrcntial operator which links 
generalized symmetries |093 

(4.4) G(J+'') =7^G(^'^ jeN, 

where s is the seed (s = 1 in most, but not all cases) and G'^-'^ is the j-th symmetry. 

A generalized symmetry, G(u), leaves 1)4. l|l invariant under the replacement 
u u + eG within order e. Hence, G must satisfy the linearized equation |093| 

(4.5) AG==F'(u)[G], 

on solutions of (|4.1() . F'(u)[G] is the Frechet derivative of F in the direction of G. 
For details about the computation of generalized symmetries, see ,GH97h. GH99j . 
A conservation law |093| . 

(4.6) Dtpix,t) + D.^J{x,t)^0, 

vahd for solutions of 1)4.111 . hnks a conserved density p{x,t) with the associated 
flux J{x,t). For details about the computation of conservation laws, see jGH97al 
IC^H97hj. 

If ()4.1|l is scaling invariant, then its conserved densities, fluxes, generalized 
symmetries, and recursion operators are also dilation invariant. One could say they 
'inherit' the scaling symmetry of the original PDE. The existence of an infinite 
number of symmetries or an infinite number of conservation laws assures complete 
integrability j093| Once the first few densities and symmetries are computed, a 
recursion operator can be constructed with the following algorithm. 



12 



DOUGLAS BALDWIN, WILLY HEREMAN, AND JACK SAYERS 



4.2. Algorithm and implementation. 

Step 1 (Determine the rank of the recursion operator). The rank of the recur- 
sion operator is determined by the difference in ranks of the generaUzed symmetries 
it hnks, 

(4.7) rank 7^y = rank cf - rank Gf\ 

where 7?, is an M x M matrix and G has AI components. 

Step 2 (Determine the form of the recursion operator). The recursion operator 
naturaUy spUts into two pieces |HG99| . 

(4.8) 7^ = 7^o + 7^l, 

where TZq is a differential operator and TZi is an integral operator. The differential 
operator, TZq, is a linear combination (with constant coefficients) of terms of type 
Dl.u^ S Z"*"), which must be of the correct rank. To standardize the form of 
TZq, propagate to the right, for example, D'^u = U2xl + u^Dr,. + uD^. 
The integral operator, TZi, is composed of the terms 

(4.9) ^^c,,G«i?-i®p(-'-)' 

i 3 

of the correct rank, where ® is the matrix outer product and p^-'^* is the covariant 
(Frechet derivative of p'--'-'). To standardize the form of 7?.i, propagate Dx to the 
left, for example, D~^UxDx — Uxl — D^^U2xI- 

Step 3 (Determine the coefficients). To determine the coefficients in the form 
of the recursion operator, we substitute TZ into the defining equation HG99, W98], 

(4.10) ^ + n'[F] + n o F'(u) - F'(u) on = o, 

where o denotes a composition of operators, TZ'[F] is the Frechet derivative of TZ in 
the direction of F, and F'(u) is the Frechet derivative with entries 



k 



4.3. Examples of scalar and matrix recursion operators. 

Example 4.2 (Kaup-Kupershmidt). Using the weights w{u) = 2,w{Dx) = 1, 
and w{Dt) = 5, we find 

rank G^^^ = 3, rank G'-^^ = 7, rank G^^^ = 9, 

(4 12) 

rank G^'^'^ = 13, rank G*^) = 15, rank G'*^' = 19. 

We guess that rank TZ — 6 and s — 2, since rank G^^^ — rank G^^-* ^ rank G^'^-' — 
rank G^^) but rank G^^) - rank G^^) = rank G^"*) - rank G^^) = 6. 
Thus, taking all -D^w^ E Z+) such that rank -D^w^ = 6 gives 

(4.13) 7^o = ciDl + C2uDl + csUxDI + au^Dl + C5U2xDI 

+ CaUUxDx + CjU^xDx + C^U^I + Cgull + CiaUU2xI + CiiU^xI- 
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Using the densities p'^^-' — u and p'^' — Su'^—Au'^, and the symmetries G*^^^ = Ux, and 
G^"^^ =F{u) — 5u'^Ux + ^UxU2x + 5uu3x+U5x from (|1.2() . we compute 

(4.14) = cuUx[D-\u2xI + 2u^I) - uxl] + cuG'-^^D~\ 

Substituting TL = TZq+TZi and G^^^ = F into H4.10|l gives us 49 linear equations 
for Ci- Solving, we find 

4c9 8c9 24c9 12c9 98c9 40c9 

. ,^ " 69 ' " " " " ~69~' " 

^ ^ _ 70c9 _ 16c9 _ 82c9 _ 26c9 , _ 2c9 _ 4c9 

~ "69"' ~ "69"' ~ "eg"'""" ~ "eg"'""'' ~ 69 ''''^ ~ "69 ' 

where C9 is arbitrary. Taking C9 = 69/4. we find the recursion operator in |W98| : 
(4.16) n = Dl + QuDl + l^xDl + Qu^Dl 

+ Y'^2xDl + SOuUxDx + Y'^3xDx + 4u^I + ^ull 

+ Y^U2xI + yU4x/ + ^UxD-\u2x + 2u^)I + G^^'>D-\ 

Example 4.3 (Hirota-Satsuma) . Only when a — \ does p.3|l have infinitely 
many densities and symmetries. The first few are 
(4.17) 

pW^^.,p(^)^3.2-2^;^ = f ■ G^^^ ^ f 5 \ ^ f a{&uux+u,x)-2vvx \ 

\ Vx J ' \ F2 J \ -[iuVx+V^x] J 

We also computed the G'^'^^ and G^^^ , but they are not shown due to length. Solving 
the weight equations 



(4.18) 



w{u) + w{Dt) = 2w{u) + 1 = w{u) + 3 = 2w{v) + 1, 
w{v) + 'w{Dt) — w{u) + w[v) + 1 = w{v) + 3, 



yields w{u) ~ w{v) — 2 and w{Dt) = 3. Based on these weights, rank p^^^ = 
2, rank p'^' = 4, and 

rank G^i' ^ ( l) , rank G^^) = f ^ 



rank G^^) = f ! V rank G^^) ' ^ 



(4.19) 

7 y \^ 9 

We would first guess that rank TZij — 2 and s = 1. If indeed the symmetries were 
linked consecutively, then 

(4 20) TZ ^ ( '^^^^ ^ ^'^^^ '^^^^ "^^^^ ^ '^^"^ ^ "^^^^ 

° \ CjDl + CsUl + Cgvl CiqDI + Ciiul + Ci2Vl 

Using H4.17|l . we have 
7^l = cnG(i)i?-i ® p(i)' = an ( ^J^^ ) i?"^ ® ( / ) = c,, ( ^f^X \ 
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Substituting 7^ = 7^o +7^l into ^J^, we find ci ~ ■ ■ ■ ~ C13 = 0. Therefore, either 
the form of TZ is incorrect or the system does not have a recursion operator. Let us 
now repeat the process taking s = 2, so rank TZij = 4. Then, 

(4.21) n = ( ) + c^,G('^D-' ® + S2iG(^)2?-i ® 

V l'^0j21 ['<-0)22 ) 

with {TZojij a hnear combination of {I?^, uD^, vD^,UxDx, VxDx,u^, uv, v'^,U2x, V2x}- 
For instance, 

(4.22) (7^o)l2 = cnDi + ci2uDl + c.^vDI + cuu^Dx 

+ Cl5VxDx + CieU^I + Cnuvl + Cisv'^l + CigU2xI + C20V2xI- 

Using lUTTIl . the first term of 7?.i in is 

7^« = 5i2G(i)7^-i ® ^ (^^- ^ ( 6ul -AvI ) 

_ ( 3uxD-^uI ~2uxD-^vI \ 

The second term of 7?.i in H4.21|l is 

Substituting the form of 71 = 7^o + 7^l = 7^o + n[^^ + Tlf' into the hnear 

system for q has a non-trivial solution. Solving the linear system, we finally obtain 



(4.23) 7^ = 

where 



7^ll 7^l2 

7?.21 7?.22 
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7^n =13^ + SuL*^ + YlUxDx + IGm^/ + 81*2:^/ - —v^l 

+ \UxD-^ul + 12MU:^i:i-i + 2u3^L)-i - ^WxD-^, 
20 . 16 ^1 16 ^ 4 ^ 8 ^ 

7^2l = -IQvxDl - I2v2xl + AvxD-^uI - Uuv^D-^ ~ AvaxD'^ 

16 8 

7^22 = -41?^ - 16uDl - 8uxDl - — u^/ - -VxD^^vI. 

o o 

A similar algorithm for the symbolic computation of recursion operators of 
systems of differential-difference equations (DDEs) is given elsewhere in these pro- 
ceedings |HSSW04j . The full implementation of these algorithms for PDEs and 
DDEs is a work in progress. 



References 

[AC91] M. J. Ablowitz, P. A. Clarkson, Solitons, Nonlinear Evolution Equations and Inverse 
Scattering, London Math. Soc. Lect. Note Ser. 149 (1991), Cambridge (U.K.): Cam- 
bridge University Press. 

[ARS80] M. J. Ablowitz and A. Ramani and H. Segur, A connection between nonlinear evolution 
equations and ordinary differential equations of P-type. I. & II, 3. Math. Phys. 21 
(1980), 715-721 and 1006-1015. 



PAINLEVE TEST, SPECIAL SOLUTIONS, AND RECURSION OPERATORS 



15 



[BGHOS] D. Baldwin, U. G6kta§, W. Hereman, L. Hong, R. S. Martino and J. C. Miller, Sym- 
bolic computation of exact solutions expressible in hyperbolic and elliptic functions for 
nonlinear PDEs, J. Symb. Comp. (2003). In Press. 
[BH03] D. Baldwin and W. Hereman, The Mathematica packages are available at 

"http:// www . mines. edu/fs_home/ whereman / 
[COO] A. R. Chowdhury, Painleve Analysis and its Applications, Monographs and Surveys in 

Pure and Applied Mathematics 105 (2000), Boca Raton: Chapman & Hall/CRC. 
[C99] R. Conte (ed.). The Painleve Property: One Century Later, CRM Series in Mathemat- 
ical Physics (1999), New York: Springer- Verlag. 
[GH97a] U. Gokta^ and W. Hereman, Symbolic computation of conserved densities for systems 
of nonlinear evolution equations, J. Symb. Comp. 24 (1997), 591-621. 

[GH97b] , Mathematica package InvariantSymmetries.m is available since 1997 from 

http://w ww.mathsource.com/cgi-bin/msitem70208-932 

[GH99] , Algorithmic computation of higher-order symmetries for nonlinear evolution 

and lattice equations. Adv. Comp. Math. 11 (1999), 55—80. 
[HG99] W. Hereman and U. Goktaf^, Integrability Tests for Nonlinear Evolution Equations. 

Computer Algebra Systems: A Practical Guide, (M. Wester, ed.), Wiley, New York, 
1999, pp. 211-232. 

[HGCM98] W. Hereman, U. Gokta§, M. D. Colagrosso and A. J. Miller, Algorithmic integrability 
tests for nonlinear differential and lattice equations, Comp. Phys. Comm. 115 (1998), 
428-446. 

[HN97] W. Hereman and A. Nuseir, Symbolic methods to construct exact solutions of nonlinear 

partial differential equations. Math. Comp. Sim. 43 (1997), 13—27. 
[HSSW04] W. Hereman, J. A. Sanders, J. Sayers and J. P. Wang, Symbolic computation of 

conserved densities, generalized symmetries, and recursion operators for nonlinear 

differential- difference equations, present proceedings. 
[HT90] W. Hereman and M. Takaoka, Solitary wave solutions of nonlinear evolution and uiave 

equations using Macsyma, J. Phys. A: Math. Gen. 23 (1990), 4805-4822. 
[HZ97] W. Hereman and W. Zhang, Symbolic software for soliton theory. Acta Appl. Math. 

39 (1995), 316-378. 

[J03] R,. S. Johnson, On solutions of the Camass-Holm equation, Proc. Roy. Soc. Lond. A 
459 (2003), 1687-1708. 

[KJH97] M. D. Kruskal and N. Joshi and R. Halburd, Analytic and Asymptotic Methods for Non- 
linear Singularity Analysis: a Review and Extensions of Tests for the Painleve Prop- 
erty, Integrability and nonlinear systems. Lecture notes in physics, (Springer- Verlag, 
Heidelberg, 1997), pp. 171-205. 

[MH96] W. Malfliet and W. Hereman, The Tanh method: I. Exact solutions of nonlinear evo- 
lution and wave equations, Phys. Scripta 54 (1996), 563-568. 

[093] P. J. Olver, Applications of Lie groups to Differential Equations, 2nd ed. Graduate 
Texts in Mathematics 107 (1993), New York: Springer- Verlag. 

[W98] J. P. Wang, Symmetries and conservation laws of evolution equations, Ph.D. thesis, 
Thomas Stieltjes Institute for Mathematics, Amsterdam, September 1998. 

[WTC83] J. Weiss and M. Tabor and G. Carnevale, The Painleve property for partial differential 
equations, J. Math. Phys. 24 (1983), 522-526. 

Department of Mathematical and Computer Sciences, Colorado School of Mines, 
Golden, CO 80401-1887, U.S.A. 

E-mail address: dbaldwinOMlnes.EDU 
E-mail address: wheremanOMlnes.EDU 



Physics Department 103-33, California Institute of Technology, Pasadena, CA 
91125, U.S.A. 

E-mail address: jackOits . caltech . edu 



